1. Field of the Invention
The invention relates to analysis of seismic trace data in petroleum exploration. In particular, the invention is a continuous method for obtaining a frequency spectrum for each time sample of a seismic trace.
2. Description of the Related Art
In the oil and gas industry, geophysical seismic data analysis and processing techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. Generally, a seismic energy source is used to generate a seismic signal that propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different acoustic impedances). The reflections are recorded in a geophysical time series by seismic detectors located at or near the surface of the earth, in a body of water, or at known depths in boreholes, and the resulting seismic data may be processed to yield information relating to the location of the subsurface reflectors and the physical properties of the subsurface formations.
Seismic or acoustic energy may be from an explosive, implosive, swept-frequency (chirp) or random source. A geophysical time series recording of the acoustic reflection and refraction wavefronts that travel from the source to a receiver is used to produce a seismic field record. Variations in the travel times of reflection and refraction events in these field records indicate the position of reflection surfaces within the earth. The analysis and correlation of events in one or more field records in seismic data processing produces an acoustic image that demonstrates subsurface structure. The acoustic images may be used to aid the search for and exploitation of valuable mineral deposits.
Accurately determining the localized spectrum of a time series or seismic trace has been a dream of geophysicists. Various techniques have been utilized in time-frequency analysis in the prior art. Traditionally, the Fast Fourier Transform (FFT) and discrete Fourier transform (DFT) have been applied. Both techniques have limited vertical resolution because the seismogram must be windowed. The spectral energy is distributed in time over the length of the window, thereby limiting resolution.
A conventional way to determine a localized spectrum is to use the method of the Short Time Fourier Transform (STFT). In STFT, tapered moving windows of the time domain signal are used to compute their Fourier spectra. The general form of this transform is:             STFT      τ        ⁡          (      f      )        =            ∫              -        ∞            ∞        ⁢                  f        ⁡                  (          t          )                    ⁢              g        ⁡                  (                      t            -            τ                    )                    ⁢              ⅇ                                            -              ⅈ                        ⁢                                                   ⁢            ω            ⁢                                                   ⁢            t                    ⁢                                                     ⁢              ⅆ        t            where f(t) is the time-domain function, g(t) is the window function, and e−iωt is the Fourier kernel. This method has been summarized by Nawab and Qauatieri (see Nawab et al., 1988. “Short-time Fourier transform” in Advanced Topics in Signal Processing. Lim, J., and Oppenheim, Al. Eds.: Prentice Hall Signal Processing Series, 289-337.) and employed in practice by Partyka et al. (see Partyka et al. 1999. “Interpretational aspects of spectral decomposition in reservoir characterization: The Leading Edge, 18, 353-360.) and Marfurt at al. (see Marfurt et al., 2001. “Narrow-band spectral analysis and thin-bed tuning: Geophysics, 66, 1274-283.) The analysis window function plays an important role in the STFT, where resolution issues are dependent to a large degree on window size. The longer the window size is in time, the better the resolution of the local spectrum in the frequency domain, but the worse the resolution in the time domain. The shorter the window size is, the better the resolution of the local spectrum in the time domain, but the worse the resolution in the frequency domain.
To solve the window function problem and improve the resolution in frequency domain, Burg (see Burg, J. P., 1968: Maximum entropy spectral analysis. Modern Spectrum Analysis, 34-48. IEEE Press.) proposed a method commonly known as the Maximum Entropy Method (MEM) and sometimes termed Autoregressive (AR) spectral analysis. MEM methods can also be found in Marple (see Marple, L., 1980, A new autoregressive spectrum analysis algorithm: IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-28, 441-454.). MEM assumes that the sampled process is an auto-regression or satisfies the Maximum Entropy Theorem. That is, the infinite series representing the process can be predicted from finite samples observed using linear prediction or the AR model. If one can determine the AR coefficients and the prediction error, then it is possible to determine the infinite series and thus its spectrum.
In one method, for an observed sequence       x    n    =            -                        ∑                      m            =            1                    M                ⁢                              a                          M              ,              m                                ⁢                      x                          n              -              m                                            +          ɛ      n      where αM,m is the AR parameter m of the Mth order AR process. εn represents the noise components. In terms of vectors,εn=XMT,nAmwhere       A    M    T    =                    [                  1          ,                                           ⁢                      a                          m              ,              1                                ,          …          ⁢                                           ,                      a                          M              ,              M                                      ]            ⁢                           ⁢              X        M        T              =          [                        x          n                ,                  x                      n            -            1                          ,        …        ⁢                                   ,                  x                      n            -            M                              ]      and T denotes the vector transpose. If each side of the above equation is pre-multiplied by complex conjugate vector X*M,n and the expected value taken, thenΦMAM=PMwhere                               Φ          M                =                ⁢                  E          ⁡                      [                                                                                X                                          M                      ,                      n                                        *                                                                                        X                                          M                      ,                      n                                        T                                                                        ]                                                  =                ⁢                              [                                                                                ϕ                    0                                                                                        ϕ                    1                                                                    ⋯                                                                      ϕ                    M                                                                                                                    ϕ                    1                    *                                                                                        ϕ                    0                                                                    ⋯                                                                      ϕ                                          M                      -                      1                                                                                                                    ⋮                                                  ⋮                                                                                                                                           ⋮                                                                                                  ϕ                    M                    *                                                                                        ϕ                                          M                      -                      1                                        *                                                                    ⋯                                                                      ϕ                    0                                                                        ]                    =                                    (                              M                +                1                            )                        ⁢                          (                              M                +                1                            )                                          is the Toeplitz autocorrelation matrix, φ1=E(xjx*j+t) is the autocorrelation function at lag time t. PM=[Pm, 0, . . . ,0]t and pM=E[εnε*n]T is the white noise power spectral density.
By definition, the white noise is uncorrelated with all xm for m<n. Thus E[εnx*n]=E[εnε*n].
It is proven that the power spectral density Sx(f) of the series is related to the input noise power spectral density by             S      x        ⁡          (      f      )        =            p      m                                        1          +                                    ∑                              m                =                1                            M                        ⁢                                          a                                  M                  ,                  m                                            ⁢                              exp                ⁡                                  (                                                            -                      ⅈ2π                                        ⁢                                                                                   ⁢                    f                    ⁢                                                                                   ⁢                    m                    ⁢                                                                                   ⁢                    Δ                    ⁢                                                                                   ⁢                    t                                    )                                                                                2      
MEM can achieve excellent frequency resolution but can be unreliable if the signal violates the assumptions of the method or if the window is too short. The major disadvantage is that it seems to be unstable, especially for less-than-expert users.
The latest approach to precluding or diminishing the problems due to windowing involves the use of the wavelet analysis. Wavelet analysis is a newly established (since the late 1980s) field in mathematics and signal processing. Like the Fourier transform, the wavelet transform also convolves through a discrete summation or continuous integration the time function (signal) with a kernel function. A method of using wavelets is found in Chakraborty et al. (see Chakraborty et al, 1995, Frequency-time decomposition of seismic data using wavelet-based methods: Geophysics, 60, 1906-1916.) as well as Xia (Xia, L., 1999, Spectra analysis of seismic data using wavelet transform: M. S. Thesis, University of Oklahoma.). In wavelet analysis, a wavelet is used as the kernel function in place of the Fourier kernel. For example, given a function f(t), its Fourier transform is:                               F          ⁡                      (            f            )                          =                              ∫                          -              ∞                        ∞                    ⁢                                    f              ⁡                              (                t                )                                      ⁢                          ⅇ                                                -                  ⅈ                                ⁢                                                                   ⁢                ω                ⁢                                                                   ⁢                t                                      ⁢                          ⅆ              t                                                          (        2        )            where ∩=2πf and e−iωt is the Fourier kernel. Its corresponding wavelet transform is:             F      ω        ⁡          (              σ        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  f        ⁡                  (          t          )                    ⁢                                    ψ                          σ              ,              τ                                ⁡                      (            t            )                          _            ⁢              ⅆ        t            where ψ is the complex conjugate of wavelet ψ. Wavelet ψ is a function that is square-integrable (ψ∈Lt(R) having zero mean and localized in both time and frequency. Therefore, ψ represents a family of wavelets that satisfies the conditions:             ψ              σ        ,        τ              ⁡          (      t      )        =            1              σ              ⁢                   ⁢    ψ    ⁢                   ⁢          (                        t          -          τ                σ            )      where τ≠0 and σ≠1. σ is referred to as the scale or the wavelet find τ is referred to as the translation parameter. Note that the wavelet is normalized such that ∥ψ∥=1.
Significantly, since wavelets are required to be localized both in time and frequency, they do not have the windowing problem of the Fourier transform and are thus well suited for localized spectral analysis. There are several different ways to cast the time series into wavelets and to compute spectral distribution therefrom. A typical method developed by Mallat et al. (see Mallat et al., 1993. “Matching pursuits with time-frequency dictionaries: IEEE Transactions on Signal Processing, 41, 3397-3415.) especially for localized spectral analysis or time-frequency analysis of seismic data is known as Matching Pursuit Decomposition (MPD).
In MPD, a family of wavelets is defined by the form             ψ              (                  σ          ,          ξ          ,          τ                )              ⁡          (      t      )        =            1              σ              ⁢                   ⁢    ψ    ⁢                   ⁢          (                        t          -          τ                σ            )        ⁢          ⅇ              ⅈ        ⁢                                   ⁢        ξ        ⁢                                   ⁢        t            where τ≠0 and σ≠0. σ a is referred to as the scale, τ as the translation parameter, and ξ as the frequency modulation. Each wavelet in the family is called a time-frequency atom. If Ω(t) is Ganssian, these atoms are called Gabor atoms. As shown by Mallat et al., Gabor atoms provide excellent time-frequency resolution. These atoms have combinations of all possible time and frequency widths and as a result constitute a redundant set. Once atoms are defined, a best match between the signal and these atoms is found by projecting the atoms onto the signal and then computing the maximum. A residue is then computed by subtracting from the original signal the product of the atoms and the cross product of the selected atom and the signal. This decomposition is continued until the energy of the residue falls below some threshold. This method has the ability to obtain a good resolution in both time and frequency for data within an intermediate frequency range like seismic data.
Wavelet transformations have been used in filtering of seismic traces in prior art. In U.S. Pat. No. 6,556,922 B2, issued to Anno, a method is described of designing and applying filters for geophysical time series data that comprises obtaining a plurality of spatially related geophysical time series and transforming the time series using a wavelet transformation. The wavelet transformation coefficients may be organized into a plurality of sub-band traces. The method includes modifying one or more transform coefficients within a plurality of the sub-band traces or within all but one of the sub-bands of traces and then inverting the sub-band traces using an inverse transform to produce a filtered version of the transformed portion of the geophysical time series. The method allows for design and analysis of non-stationary filters and filter parameters in untransformed data space. Non-stationary signals may be filtered in all or in windowed portions of the geophysical time series data.
Interference between wavelets is a key problem with wavelet transform methods. It is possible to obtain a very good spectrum for a whole time series using the wavelet transform, without concerning oneself with interference effects. But for local spectral computation, interference effects cannot be overlooked. There is a need for a localized spectral analysis method in seismic data that gives an instantaneous spectrum of the seismic signal at each instant of time. The present invention fulfills this need.